(N 

o 

(N 



mm ipPlEg^ Published by Institute of Physics Publishing for SISSA/ISAS 



Received: 
Accepted: 



Direct recovery of density fluctuation spectra from 
>>: tomographic shear spectra 



Pino Torinese (Torino), Italy 



Keywords: cosmology: theory, dark matter, gravitation; methods: numerical, N-body 



simulations. 



M 
O 



Marino Mezzetti 1,2 , Silvio A. Bonometto 1 ' 23 , Luciano Casarini 4 , Giuseppe Murante 1,5 

f) • 1 - Department of Physics, Astronomy Unit, Trieste University, Via Tiepolo 11, 

_; ! / 34143 Trieste, Italy 

2 - I.N.A.F. - Astronomical Observatory of Trieste, Via Tiepolo 11, I 34143 Trieste, 
' Italy 

3 - I. N.F.N. - Sezione di Trieste, Via Valerio, 2 I 34127 Trieste, Italy 

4 - Departamento de Fisica, UFES, Avenida Fernando Ferrari 514, Vitoria, Espirito 
Santo, Brasil 

5 - I.N.A.F. - Astronomical Observatory of Torino, Strada Osservatorio 20, 1 10025 

V ' ' Pun n T^nrinpQP iT^nTinn) Ttnhi 



ISO 



o 



> 

t> 

Abstract: Forthcoming experiments will enable us to determine high precision tomo- 
graphic shear spectra. Matter density fluctuation spectra, at various z, should then be 

£X recovered from them, in order to constrain the model and determine the DE state equa- 

tion. Available analytical expressions, however, do the opposite, enabling us to derive shear 
spectra from fluctuation spectra. Here we find the inverse expression, yielding density flue- 

S^ tuation spectra from observational tomographic shear spectra. The procedure involves 

SVD techniques for matrix inversion. We show in detail how the approach works and 
provide a few examples. 
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1. Introduction 

Dark Energy is the main finding and puzzle in today's cosmology. Its contribution to the 
cosmic budget is directly constrained by Cosmic Microwave Background spectra, but only 
measures of the matter density field p(x, z), at low z, can provide clues on its state equation 
w(z). It is then important that weak lensing data, directly sensitive to the whole matter 
distribution, can be translated into information on the density fluctuation spectrum 



P(k,z) = {\6(k,zf 
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and, namely, on its redshift (z) dependence. Here S(k, z) is the Fourier transform of the 
matter fluctuation field e(x, z) = p(x, z)jp—\. The final target would be recovering P(k, z) 
from tomographic shear data. 

As a matter of fact, data can be used to obtain the tomographic shear spectra Cij(£) - 
defined below - and a known relation yields Cij(£) from P(k, z). In this paper we therefore 
aim at inverting such a relation, so obtaining P(k,z) from Cij(£). 

The problem somehow reminds the inversion of the Limber equation (see, e.g., Q), 
aiming to obtain the angular 3-point correlation function £(r) (making no ansatz on its 



o 
o 
o 



According to Huterer & Takada |1C], data may become precise enough to enable us 
to appreciate cosmological parameter variations causing a shift 0(1%) in P(k,z). This 
conclusion was attained by parameterizing deviations in P(k, z) and testing their conse- 
quences on the angular shear spectra, trying also to take into account all possible sources 
of systematic bias, but assuming that the deviations in the shear spectra are solely due to 
the shift in the density power spectrum P(k,z). 

In our case, while the power spectrum and its evolution measure the linear and non- 
linear growth factors Q(k,z), the relation between P(k,z) and the shear spectra depends 
on a kernel involving a number of astrophysical assumptions (concerning, e.g., the galaxy 
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form) from the 2-point spatial function w(6). There are quite a few differences, however. In 
particular, 3-D galaxy catalogs made the Limber equation obsolete. On the contrary, even 
having spectroscopic redshifts for all lensed galaxies, no 3-D shear spectrum is recovered, 
although the 2-D shear spectrum would then be known with great accuracy. Furthermore, 
in order to recover fluctuation spectra, shear spectra at various redshifts are to be simulta- 
neously used; on the contrary, the inversion of the Limber equation, even in the relativistic 
regime [||, is independently effective at various redshifts. 

Many authors debated the use of tomographic shear spectra to constrain the cosmo- 
logical model, with different procedures y, [|] (see also ||). Also the Dark Energy task 
force (DEFT: ||) devoted much attention to this approach. The basic pattern essentially 
amounts to comparing observational Ctj(£) data with the theoretical Cij(l) obtainable for 
models belonging to an assigned parameter space. This Bayesian procedure has been used 
in quite a few former cases and is certainly effective, while the technique discussed here 
is not yet mature enough to compete with it. However, we see significant possibilities to 
upgrade it, to at least obtain a complementary tool. 

At present, significant shear data are already available. Cosmic shear measurements 
were obtained by using large area ground surveys (see, e.g., 0) or narrower area space data, 
characterized by high quality imaging (see, e.g., [§). However, observational campaigns 
to perform ultimate systematic mapping of tomographic cosmic shear are the basic aim of 
future missions. 

In particular, the Euclid project, a recently approved ESA mission, is devised to observe 
about half extra-galactic sky (~ 15,000 deg 2 ) from space, at a diffraction limited spatial 
resolution which would be impossible from ground ||. Euclid will also obtain medium 
resolution (R 400) spectra of 1/3 of all galaxies brighter than 22 mag, in a wavelength 
range unreachable from ground for faint galaxies above z=l. By measuring the correlations 
in the shapes of ~ 1.5 billion galaxies (> 2 orders of magnitude more than all galaxies 
in today's samples), Euclid will map weak gravitational lensing with extreme accuracy, 
yielding dj(£) with a precision 0(1 %) up to £ ~ 5000. At larger £ (up to ~ 30, 000) the 
precision gradually worsens. 

At low £, cosmic shear can be derived from linear fluctuations. Already at £ ~ 100, 
however, non linear contributions exceed ~ 1 % while, above £ ~ 500, neglecting non linear 
structures is clearly misleading. Baryonic physics causes spectral shifts already at £ ~ 300; 
they become more and more relevant at larger £ and, above £ ~ 2000, shear spectra are 
unpredictable if we neglect it. 
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number distribution as a function of z or the relation between photometric and physical red- 
shift) and a specific cosmological assumption, the time dependence of the scale factor a (or 
the redshift z). 

This is why a parameter shift affecting the growth factor Q(k, z) seldom leaves the 
kernel unaffected. It may well be that the kernel variations add up to G(k, z) variations, 
so strenthening the Huterer-Takada effect. The opposite case is however also possible. If 
mutual cancellations occur, when some specific parameter shift is considered, such specific 
parameter cannot be fixed with the claimed accuracy. This is not just a theoretical caveat, 
as a number of examples can be given. 

We stress this point also because a similar difficulty affects the inversion procedure 
described below. When supposing that the shear spectra are measured, we shall aim at 
reconstructing P{k,z) - and thence the growth factor Q(k,z) - from them, by following 
a direct analytical pattern. Clearly, two distinct time dependences arise from the model 
choice: (i) background equations rule the dependence of the scale factor a on time and, 
therefore, space time geometry; (ii) fluctuation dynamics rules the time dependence of the 
growth factor Q. To achieve our dynamical aim, the procedure we describe will assume 
that the geometrical kernel is assigned. 

In the discussion Section, we shall however return to this point. As a matter of fact, 
geometry depends only on a part of the parameters defining the model, while dynamical 
data open a window, e.g., on the separate contributions of baryons and Dark Matter to 
cosmic matter and on further parameters unrelated to the background metric, but critical 
to define primeval fluctuations and their later evolution. Furthermore, in the discussion 
Section we shall conjecture that this apparent difficulty might be turned into a tool for 
model discrimination. 

The plan of the paper is as follows: In the next Section we shall discuss the expression 
yielding Cij(l) from P(k,z); this will enable us to outline how such an expression can be 
formally inverted. In Section 3 we shall enter into technical details concerning the inversion 
procedure. In particular we shall introduce the SVD technique, essential for dealing with 
(nearly-)singular matrices. The technique will enable us to provide a concrete solution to 
the inversion problem, attaining a precision 0(1 : 1000), at least, that we shall illustrate 
by using Halofit fluctuation spectra. In the same Section, however, a number of difficulties 
will also be outlined. In Section 4 we shall go beyond Halofit, applying the technique 
to hydrodynamical simulation outputs; we shall also show an approach enabling us to 
overcome some of the difficulties previously outlined. Finally, Section 5 is devoted to a 
discussion of the results obtained and the perspectives opened. 

2. From fluctuation to shear spectra and viceversa 

The convergence weak lensing power spectra are linear functionals of matter power spectra 
at various z, suitably convoluted with the lensing properties of space, mostly due to the 
matter distribution in it, and the background galaxy distribution. We set the galaxies, 
whose images can be distorted by gravitational lensing, into n bins at increasing depth, 
labeled by i,j = l,...,n . Their distributions will limit the functions Wi(u), gauging the 
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effects of the lensing systems, in the expressions 
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Cij(£) = tf 4 / du Wi(u)Wj(u) P(£/u,u) 
Jo 



(2.1) 



yielding the tomographic shear spectra Cij{£) (T^, OL O]. Here ro is the conformal age of 
the Universe, t = tq — u being the conformal time in the FRW metric 



ds 2 = a 2 (r)(dT 2 -dX 2 ) 



(2.2) 



so that d\ 2 is the co-moving 3-space metric, that we assume to be flat; P(k,u) is the 
fluctuation spectrum at the conformal time set by u; Hq is the Hubble parameter. The 
window functions W% are then defined in the next subsection. 

2.1 Window functions 

In the literature, n = 1, 3 or 5 bins were considered. Data now available could not be 
analysed with > 5 bins. With the use of the ordinary bayesian procedure, in fact, shot 
noise in data allows us no improvement in parameter determination already when going 
from 3 to 5 bins (see, e.g., [0]). As a matter of fact, shot noise adds to the Cij(£) spectra in 
a way oc n/N g (N g : total number of lensed galaxies observed) and, above n = 3, its effects 
confuse the signal dependence on the band (i,j), namely because of the approximation 
implicit in using photometric redshifts. However, if N g increases by a factor 100, as in 
future Euclid data, even with 10 bins shot noise can be ignored. Here we shall however 
keep to a 5-bin case, to prevent numerical complications in matrix algebra. 

The bin limits z% are conveniently selected so to have the same number of galaxies per 
bin; we shall also assume, as usual, that the distribution of the galaxy number in redshift 
and solid angle reads 
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and A = 2, B = 1.5, so that C = 1.5/zo (with zq = z m /lA12 obtained from the median 
redshift z m = 0.9 ). 

This distribution is then considered within the limits of the redshift bins, taking how- 
ever into account that only photometric redshift values are given. The discrepancies be- 
tween them and the actual galaxy redshift define the filters 
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Figure 1: Distribution of galaxy redshift values and actual shape of (5) bins, if denned by using 
photometric redshift. In the inner frame the z-u relation is shown (u in Mpc). The plots are for a 
spatially flat model with £7 m = 0.24 and Hq = 73 km/s/Mpc. 
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Figure 2: Wi(z) functions, yielding the effective distribution of lensing systems, in the 5-bin case. 

with a(z) = 0.05 (1 + z) coherently with Euclid expectations || (see also 1 1]]) and set 

Di{z) = n{z)Iii{z) (2.6) 

yielding the distributions 
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Figure 3: To obtain shear spectra, the fluctuation spectra P(k,u) are to be integrated. The lines 

are examples of the domains of integration, for logarithmically equispaced i = 4000, 8000, , 

32000 , on the plane spun by k and u = t — r. 



Si(z) 



Di(z) 



(2.7) 



f °° D t (z')dz' 

as a function of z of the actual setting of lensed objects. Figure [j] exhibiting the resulting 
redshift bins, shows that increasing their number causes extensive overlaps, which risk 
polluting the analysis when the number of galaxies per bin is too small. 

Within Figure [l] we also show the conversion between z and u (in Mpc), for a specific 
model consistent with WMAP-7 CMB results |0§, with fi m = 0.24 and H = 73 km/s/Mpc 
(matter density and Hubble parameter). This model will be used throughout this work, 
in order to exploit the results of wide simulations with a large dynamical range, including 
also baryon physics, available to us. 

The images of any galaxy belonging to a bin can be lensed by systems at lower z. From 
the distributions SAz) we therefore derive the functions 
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u(z') 
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The factor in square brackets must be set to zero if negative: the lens is closer than the 
lensed galaxy. Such functions then yield the window functions 

3. 



Wi{z) 



-n m F t {z)(i + z) 



(2.9) 



to be used in eq. (|2.1|), In Figure || we show the Wi(z) profiles in the 5-bin case. 

2.2 Fluctuation spectra 

The dependence of Cy on 



aperture $ 



can be roughly interpreted as a dependence on an angular 
2ir/£ which subtends linear scales A = 9u, increasing with z. In turn, the 
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rough correspondence k ~ 2tt/X also holds. Altogether, a given £ corresponds to decreasing 
k values as u (or z) increases. 

Accordingly, in eqs. (|2.1| ), the fluctuation spectrum is taken at decreasing k values, as 
u increases. In principle, for u — > 0, P(k) should be evaluated at k — > oo where, however, 
it vanishes. 



If we consider the logfc-logu plane, the integration in eq. (2A) is carried along tilted 
straight lines such as those shown in Figure |3|, each line corresponding to a given £. 

The spectra P(k,z) therefore yield functions Pg(u), where £ fixes a line on the logfc- 
logu plane and u is used as an abscissa for such a line. More explicitly, it will be Pt{u) = 
P[£/u,z(u)]. 

2.3 From fluctuation to shear spectra 

The integration interval in eq. (2.1) is apparently finite. When u reaches the conformal 



age of the Universe To, however, z approaches oo. Figure shows that all Wj vanish well 
before so, and this sets an effective upper limit to the integration. 

Aiming at 6— digit precision, a numerical integration performed by summing on a large 
number of equispaced points requires ~ 5000 points up to u ~ 6000 (z ~ 3). By itself, 
however, a large number of points does not guarantee a safe result; what matters is that 
Pi{u) has a fair value at each u considered. The point is that Pe(u), for each £, might result 
from interpolating along a set of points u r and the key to obtain fair Cij(£) values is that 
the number of u r is adequate. 
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Two cases are considered here: we first used Halofit [16], enabling us to evaluate 
approximated spectra at any k = l/u for any u(z); then, we used more precise spectra 
worked out from simulations, although at a given set of redshifts. 

In principle, in the former case, the Halofit package can be directly questioned for 
each k and u. Vice-versa, Halofit can also be exploited to test how many k s and u r 
values are needed to obtain interpolated results equivalent to those resulting from "direct 
questioning" . 

A point one appreciates soon is that interpolating along the tilted lines of Figure || 
is not so effective as interpolating among spectra at constant k along u (or z). More 
quantitatively, Halofit was used to test that using the results of a large simulation, whose 
fluctuation spectra are known at the redshifts 

l + z r = 10 r/20 (r = l, ....,19) (2.10) 

is adequate to achieve the required precision. This is true if we interpolate at constant k. 
In contrast, if we deduce first P^(u r ) values and interpolate then among them, fluctuation 
spectra at more z r (approximately 3 times as many) are needed. 

The cosmology for which we consider Halofit spectra is the same used to run the 
simulation in the next Section. More specifically, we take a flat ACDM model with f2 m = 
0.24, Slj, = 4.13 x 10~ 2 , h = 0.73, n s = 0.96, (density parameters of total matter and 
baryons, Hubble parameter, primordial spectral index, respectively) and normalized so 
that the m.s.a. of density fluctuations, at 8/i _1 Mpc, as = 0.8 at z = . 



For our purposes it is necessary to make also use of Gaussian integration, i.e. to project 
the integrand function f{x) onto polynomials Tt a (x), orthogonal with an assigned weight 
function R(u), finding its components f a . The integrals H a of each polynomial are then 
known and Yua=l fcJ^-a is a reliable integral of f(x), if N is large enough. As is known, 
this technique can be translated into a practical and simple procedure, so that integration 
is reduced to a weighted sum of values taken by the integrand function f(x) in a suitable 
set of points x a . 

More in detail, using monic polynomials, we have that 

dx R(x)ir a (x)Trp(x) =N5 aj3 , (2.11) 

o 

with a known normalization J\f . Monic polynomials are obtained from a suitable recurrence 
relation assuming that the coefficient of the leading term, for each a, is unity. If we 
then truncate the sum to N terms, the N zero's of ttn( x ) are the points x a , while the 
corresponding weights are 
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f£°dx Rjx^^jx) 

TTN-l(x a )lT' N (x a ) 

tt'(x) being the ordinary derivative of ir(x). Then, 



(2.12) 



dxf(x) = J2 W af{Xa) (2-13) 



r=l r=l 



A = ij, c A = C ij /H* (2.15) 



with the correspondence law 



while 



i,j 1,1 ... 1,5 2,2 ... 2,5 3,3 ... 5,5 
A 1 ... 5 6 ... 9 10 ... 15 



S A , Xr = Wiix^Wjix^/Rixr) , (2.16) 



and, for R(x) oc e~ x , this technique is dubbed Gauss-Laguerre integration, as ir a (x) = 
L a (x), the Laguerre polynomials. 

In the case of eq. (^J), where integration is cut off by the exponential-like decay 
of the Wi functions, this approach can be applied by assuming that u = 4>(x) or, more 
specifically, x = (u/u) 13 and suitably selecting then u and (3, 

We shall show the degree of approximation allowed by such a technique after discussing, 
in the next subsection, why it is needed and which are the limitations to N. 

2.4 Formal inversion 

We then rewrite eq. ( |2,1D as follows: 

N N 

C A{€) =Y / WrS A ,x r Px r (t) = Y, M ^Pr(i) ■ (2-14) 

Here we have set 



p r (£)=p Xr (£) = P t [<l>(x r )] = P[£/<f>(x r ),(f>(x r )} . (2.17) 

If, in eq. (2.14), we take N = 15, Mat are square matrices and, provided that they are 



not singular, the inverse equation 

Px M = ^2(M);lc A (£) (2.18) 

A 

also holds. Accordingly, we shall be able to recover the spectrum P(k, z) for any k = 
£/4>(x r ), at the redshift values z[4>(x r )]- 

Notice that the inversion procedure acts on each £ value separately. Any z and/or k 
value can be attained, in principle, just by suitably choosing the u and j3 parameter, to 
obtain a suitable <j)(x) function. In principle, the choice can depend on £. 

M 



All this makes it clear that Gauss-Laguerre summations ( 2.14 ), in this case, must be 
limited to N = 15. This is a consequence of using 5 bins. With 3 bins, A^ = 6 at most, a 
value inadequate to yield any reliable integration. On the contrary if, e.g., we take 7 (10) 
bins, we have A^ = 28 (55). When increasing the number of bins, we therefore approach 
an increasingly satisfactory situation, as many terms can be set into the summation ( 2.14| ) 
and, accordingly, many more linear equations can be used. 
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In order to divide the galaxy sample into many redshift bins, however, we must have 
either reliable redshift values for most lensed galaxies, or a sample including very many 
galaxies. The first option probably requires one to do better than using photometric 
redshifts; the latter option is the one pursued by Euclid. In this paper, however, we shall 
consider only the 5-bin case. 

3. Operational problems 

When trying to exploit this formal inversion we find two kinds of difficulties: (i) the 
performance of the integration procedure; (ii) a quasi-singular behavior of the matrix Mat- 
Both of them can be, at least partially, overcome and the results we give here aim to 
show that the procedure is effective. 

3.1 Integrations 

For instance, in order to avoid a Mat singular behavior, a possible option is to reduce its 
dimension. Before discussing this in more detail, we discuss the results of a Gauss-Laguerre 
integration when reduced to 12 points. 

In Figure || we show £ 1 - 2 ca{£) (A = 1, ... 15) obtained by using Halofit spectral expres- 
sions. £ values from 10 to 5000, at intervals of 5, are taken. Three integration techniques 
are compared: (i) The black curve is the benchmark obtained by performing a Riemann 
integration with 10000 points between u = and 6000 Mpc. (ii) The magenta curve is 
the Gauss-Laguerre integral with u = 1200 and (5 = 2.9. (hi) The green curve is also the 
Gauss-Laguerre integral with u = 600 and /3 = 1.9 . 

This Figure allows us to draw some immediate conclusions: (a) The performance of 
a Gauss-Laguerre integration exhibits a strong dependence on the parameter choice, (b) 
Results are however better for larger i, j. 




LogW 

Figure 4: Tomographic shear spectra obtained from Halofit fluctuation spectra. We plot CV, x I 1 - 2 , 
so to reduce the ordinate range and stress spectral discrepancies. Black spectra from Riemann 
integration (10000 points). Green and magenta spectra from the Gauss-Laguerre integration with 
different choices of u and j3 (see text). Values I = 10, 15, 20, .... , 5000 are plotted. Spectra are 
larger for greater i, j. Solid lines are Cu (i — 1,...,5). Dashed lines are CV,- (j > i). Notice 
the gradual decrease of discrepancies towards greater i, j, due to a wider z -range contributing to 
integration, i.e., to a richer data system. 
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The reason for point (b) is soon evident: The redshifts z(u r ) = z[(j>(x r )], shown in 
Table I, are distributed between and ~ 1.5. All of them yield a substantial contribution 
to the sum for large A. On the contrary, owing to the fast cut off of Wi, for small i values, 
at low A the sum risks including only a few terms yielding a real contribution. This makes 
it also clear that using u values corresponding to z >~ 1.5 is vain. The contribution to the 
integrals coming from larger z are to be approximated by the assumed exponential decay 
of the integrand function and the shape of the orthogonal polynomials selected. 

It is then also clear why the (iii) integration procedure yields better results than the 
(ii) one: In the former case, the u r = 4){x r ) values are more densely accumulated at lower z, 
thus allowing better integration for low A. Apparently, this causes no detriment to high-^4 
results. 
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Table I 

z r values 



(ii) case 












.0466 


.1152 


.1886 


.2664 


.3499 


.4402 


.5386 


.6473 


.7700 


.9122 


1.0850 


1.3162 


(iii) case 












.1234 


.2338 


.3338 


.4304 


.5263 


.6237 


.7245 


.8307 


.9450 


1.0714 


1.2171 


1.4004 



Of course, small shifts of j3 or u cause no substantial difference. However, although 
exploring different options with much care, we cannot exclude that better u, (3 choices exist. 
In particular, if one uses a different measure R(x) on the integration interval, different 
orthogonal polynomials follow. For instance, one could use R(x) oc e~ x yielding 7r a (x) = 
H a (x), the Hermite polynomials; or some other R(x) yielding non-tabulated polynomials. 
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3.2 The SVD technique 

Let us then consider the inversion procedure. The problem here is that the matrix 

Mat = w r Wi(x r )Wj{x r )/R(x r ) (3.1) 

tends to be singular. There is a specific analytical reason for that: when u r = x r u 
exceeds ~ 0.7 any matrix element containing W\ tends to vanish, as is evident from Figure 
|2[ A similar feature is caused by any element containing a generic W% with i ^ 5, when u r 
yields a redshift z ~ 1.5 . But, even keeping all u r below 1.5, the ratio between largest and 
smallest diagonal elements tends to be too large, even in double precision. 

Before further discussing this point, let us introduce a specific technique, allowing us 
to gauge the degree of singularity of a matrix, dubbed SVD (singular value decomposition). 

It is based on a theorem of linear algebra, stating that any real N r ® N c matrix Ai, 
with N r > N c , can be decomposed into a rows x columns product 

M = U x | diag(si) \ x V T . (3.2) 

Here V is the transpose matrix of a matrix V which, as well as U, is orthonormal, while 
s is diagonal. Apart from multiplicative factors, the decomposition is unique. If A^ = N c , 
the inverse of Ai , in general, reads 

M' 1 = V x | diag(l/si) \ x U T , (3.3) 

and the technique is also a valid numerical way to invert large matrices. 

The degree of singularity, however, can be inspected by just considering the Sj com- 
ponents. One (or more) vanishing element(s) cause the matrix to be singular. Even if it is 
not so, however, and the ratio between the greatest and smallest Si exceeds ~ 10 6 (10 12 ), 
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Figure 5: Recovery of P(k,z r ) spectra from dj(£); the values z r are yielded by x r values. Here 
the 4 lowest z values are shown. The ratio of input to recovered values never differs from unity 
more than ~ 10 -4 . 

there is no hope to invert Ai in single (double) precision. Anyhow, if a level of precision 
0(1 : 10 6 ) is to be kept, one must use double precision keeping the highest Si/sj ratio 
within ~ 10 6 . 

A fair discussion of the technique can be found in Numerical Recipes [17| or in Matrix 
Computations [|18|| , where is also discussed how this technique can be applied when iV r > N c , 
as well as what to do to find tentative solutions when the highest Si/sj ratio is too large, 
and even when some Sj vanishes. 

Here we shall not debate this approach any further. We shall just report that it allowed 
us to test a wide number of options, by selecting those exhibiting a lower level of singularity. 
Although our inspection was systematic and detailed, we cannot exclude that even more 
efficient solutions can be found. Here we wish to outline that solutions are indeed available, 
and then discussing which problems remain. 
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3.3 Inversion 

Making use of 12 Gaussian points to integrate, eq. (15) can then be considered as a linear 
system of 15 equations with 12 unknowns p r (£) (r = 1, .., 12). The very SVD technique is 
built to deal with such cases, and the very redundancy of the system is a reason for keeping 
the top Si/sj value within ~ 10 6 . 

If we try to recover the P(k,z) spectrum from the spectra Cij(£), we then have full 
success. This is shown in Figures || and y. The degree of precision is better for higher 
redshift values: we pass from ~ 0.2% at low z, to a precision 0(1 : 100000) for z = 1.316, 
the highest redshift considered. 
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Figure 6: As the previous Figure, for z r with r 
smaller, keeping well within 1:10 6 for r > 8. 
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6, 8, 10, 12. Discrepancies, for larger z r are even 



A point to outline is also the shift of the k interval explored when z increases. The 
excellent results for the highest redshift, e.g., are concentrated within an interval where 
non linearity is just approached, also because the non-linear fc-range shifts to the right at 
higher redshift. With Halofit spectra, however, one can hardly do better. 

We tested the SVD algorithm also using 12 equations, excluding either the 3 top-^4 
equation, or the equations for i = 1, j = 1,2,3. The inversion succeeds, discrepancies 
keep an acceptable level, but are indeed larger (and somewhere notably larger) than those 
obtained with the full set of equations. In the former case, the worst output/input ratio 
is at intermediate redshifts, where errors increase up to a factor 20 . In the latter case, 
the worst results concern highest redshift values, where the errors have a really significant 
increase, by a factor up to ~ 4 x 10 4 , although however keeping within 0(5 %). 

Figures |5| and |6| are one of the main results of this work. To operate the inversion, 
we made use of the very dj (£) obtained by using a Gaussian integration procedure (black 
curves in Figure |||). In turn, Gaussian-Laguerre integrals were based on spectra obtained 
from interpolating Halofit spectra at the redshifts 1 + z r = 10 r ' 20 (r = 0, 1, ..., 19), in order 
to retain homogeneity with the results obtainable from simulations (see next Section). 
Interpolation was performed at constant k values. 

The next point we whall inspect here is the "stability" of the inversion procedure. 
We do so by applying the inversion algorithm to the angular spectrum worked out by 
performing a Riemann integration. 

This operation does not yield an immediate success. The discrepancies between Cij{£) 
worked out with 12 Gaussian points or 10000 Riemann points is fairly large in some points, 
up to O(10%) and in Figure ^ we show the results of a brute inversion for large z r . The 
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Figure 7: Results of the inversion algorithm based on 12 point integration, if applied to "ex- 
act" spectra. Although the power spectrum is better approximated at high redshift, the result is 
unsatisfactory, namely because the z dependence of the power spectrum is not reproduced. 

results of this procedure are better at large redshifts, where the discrepancy between Cij(£) 
obtained through Riemann and Gauss integrations is smaller. 

It is however clear that the procedure is to be tested when using more than 12 Gaussian 
points, so reducing the discrepancy between Gaussian and Riemann integrations well below 
the expected noise level. Rather than the Gaussian-Riemann discrepancies, it will then 
become important to test the impact of random noise. 

However, still with 12 Gaussian points, there is a possible improvement of the results, 
that we shall debate after discussing the replacement of Halofit spectra with simulation 
spectra. 
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4. Beyond Halofit 

In this Section we will report the results of the same procedures, when applied to spectra 
obtained from large hydrodynamical simulations. This enables us to consider larger £ values 
and shear spectra exploring scales well inside galaxy clusters, which cannot be approached 
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if baryon physics is disregarded. We also recall that several authors ]l9|, |j| discussed 
the limitations in the Halofit reconstruction of non-linear power spectra, whose use was 
dangerously extended to models including DE with an equation of state different from 
w = — 1. 

4.1 Simulation 

The simulation used here follows the development of structures within a periodic box of 
comoving side L = 410/i~ 1 Mpc, where (2x)1024 3 particles are set. The cosmology is 
a spatially flat Gaussian ACDM model with tt m = 0.24, tt b = 4.13 x 1CT 2 , h = 0.73, 
n s = 0.96. Fluctuations are normalized so that the m.s.a. of density fluctuations, on the 
scale of 8 fo -1 Mpc, as = 0.8 at z = . The two populations of 1024 3 particles, therefore, 
have masses m c ~ 1.89 x 10 9 h~ l M Q and m b ~ 3.93 x 10 8 /i~ 1 M . 

The simulation was carried out by using the TreePM-SPH GADGET-3 code (SPH: 
smoothed-particle hydrodynamics), an improved version of the GADGET-2 code (Springel 
2005). Initial Zeldovich displacements were generated at Z{ n = 41. Gravitational forces 
were computed using a Plummer-equivalent softening which is fixed to a physical scale 
epi = 7.5/i -1 kpc from z = to z = 2, being then constant in comoving units at higher 
redshifts. 

As far as baryon physics is concerned, radiative cooling was computed for non-vanishing 
metallicity according to Sutherland & Dopita J2(|, also including heating/cooling from a 
spatially uniform and evolving UV background. Gas particles above a given threshold 
density are treated as multi-phase, so as to provide a subresolution description of the inter- 



branch (AGB) stars, as described in ref. [22]. Stars of different mass, distributed according 
to a Salpeter IMF, release metals over the time-scale determined by the corresponding 
mass-dependent life-times. Kinetic feedback is implemented by mimicking galactic ejecta 
powered by SN explosions. In these runs, galactic winds have a mass upload proportional 
to the local star-formation rate. The wind velocity is then v w = 500km/s; this corresponds 
to assuming about unit efficiency for the conversion of energy released by SN-II into kinetic 



energy for a Salpeter IMF. More detail on this simulation are provided in [23]. 

Simulations including baryon physics, aimed at evaluating shear spectra, were per- 
formed by various authors (see, e.g., [JE3J ||]). Van Daalen et al. [24] included also AGN 



feedback in their simulations. The simulation used here, however, was performed in a box 
and with a dynamical range large enough to enable us to evaluate fluctuation spectra for 
the whole range needed to compute Cij (£) up to £ ~ 40, 000-50,000 . 

Simulations including baryon physics, aimed at evaluating shear spectra, were per- 



formed by various authors (see, e.g., [[ll], y]). Van Daalen et al. [24] included also AGN 



feedback in their simulations. The simulation used here, however, owns two suitable fea- 
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stellar medium, according to the model described in ref. [21]. Within each multi-phase gas 
particle, a cold and a hot-phase coexist in pressure equilibrium, with the cold phase provid- 
ing the reservoir for star formation. Conversion of collisional gas particles into collisionless 
star particles proceeds in a stochastic way, with gas particles spawning a maximum of two 
generations of star particles. The simulation also includes a description of metal production 
from chemical enrichment contributed by SN-II and SN-Ia supernovae and asymptotic giant 
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tures: (i) it was performed in a box large enough (410 /i _1 Mpc) to allow us a direct connec- 
tion between simulation and linear spectra; (ii) a dynamical range large enough was used, 
enabling us to evaluate Cij(£) up to £ ~ 40,000-50,000 , by exploiting fluctuation spectra 
unaffected by numerical noise. 

4.2 Fluctuation spectra 



Power spectra are computed at the redshift values ( 2.10| ), by using the algorithm PMpow- 



Figure || and 10, however, show that the relation (^4|) is efficiently invertible, in prin- 
ciple, also for those £ which deal with the inner cluster structures. 

Let us finally return to using the inverse of the matrix Mav, obtained with 12 Gaus- 
sian points, to invert the "exact" C«(£) spectra. In Section 3 we already noted that the 
question will bear a completely different aspect when we are allowed to use more than 5 
bins. However, even with 5 bins, a substantial improvement is attained through a suitable 
renormalization procedure, essentially requiring the same information needed to obtain the 
window functions W{(z), i.e. the background "geometrical" model parameters. 

If we add to these parameters a value of the spectral index n s we easily evaluate 
the linear fluctuation spectra p( c,hn >(k,z) of a pure-CDM model, with the same Q m and 
#o- (Let us however soon outline that the choice of n s is almost arbitrary. In 
particular we do not need to use the specific value used in the simulations.) A 
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erM included in the PM package, courtesy of A. Klypin. 

Through a CiC procedure the algorithm assigns the density field on a uniform Cartesian 
grid starting from the particle distribution. It then calculates the spectrum through a 
FFT on an 3 [n = 2 4 A) grid; i.e. n = 1024 x 2 4 = 16384 . These large n are obtainable 
by considering a A 3 grid in a box of side L/2 4 , where simulation particles are inset, in 
points with coordinates Xij = x% — isL/2 4 (i = 1,2,3), and an integer v selected so that 
< Xij < L/2 4 . This technique provides spectra down to wavelengths slightly above the 
gravitational softening scale, being only limited by the numerical noise due to the grid 
used to set the initial conditions. The limitation would be however identical if large-n 
grids could be directly applied to the original box of side L. 

The spectra shown in Figure || are obtained by merging the linear spectrum of the 
model with the simulation spectrum. No recourse to approximate spectral expressions, like 
Halofit, is needed to interconnect the two spectral parts. More details on the techniques 
to merge the two spectra can be found in |23| (see also [^6| ) . 

4.3 Spectra of cosmic shear 
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The recovery of fluctuation spectra works on each £ value separately: spectra at different 
k values for each redshift z are obtained by relating the values taken by the 15 Cy at that 



£. The curves in Figures ^ and 10, accordingly, arise from the merging of contributions 
coming from the large set of £ values available. 

Formally, all Cy contribute to recover P(k, z) at any z. However, the physical reason 
why the recovered P(k, z) is slightly less precise at small z resides in the fact that the 
Cij{£) are almost-vanishing for the significant £ values. This is also the reason why the 
inversion procedure is more efficient when more Gaussian points fall at low z. 
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Figure 8: Tomographic shear spectra for i,j — 1, ....,5. Solid and dashed lines as in Figure 4. 
Gauss-Laguerre integration performed with case (ii) points. Spectra obtained from Halofit (red) 
are overlapped to those obtained from the hydro simulation (black). The latter extend well in the 
non-linear region. Discrepancies between the red and black curves are significant in the critical 
region of non-linearity onset and do not reduce towards great i, j . 

Riemann integration then yields the angular spectra C^' (£) from the fluctuation spectra 
p( c,lm '(k,z). In turn, the angular spectra can be tentatively inverted by using the very 
SVD "Gaussian matrix" worked out for the full model. The final output are linear spectra 
p( c > lin > (k, z) to be compared with the input linear spectra p( c,hn '(k, z), so building the 
normalization factors 

p( c > lin \k,z) 



N{k,z) 



p(c,lin)(k,z) 



(4.1) 



We shall use these to try to remove the bias, tentatively assuming that it is the same for 
the full model and the CDM-linear case. 

Let P(k, z) be the spectra obtainable by interpolating the spectra of the simulation at 
the redshifts ( p. 10 ). We obtain from them the shear spectra with either 12-point Gaussian 
or Riemann procedures. Then, we calculate the inverted Riemann spectrum P(k,z s ), by 
applying the inverse "Gaussian matrix" , as above. 

These steps give us results not too different from those shown in Figure [5]. We then 
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Figure 9: Recovery of fluctuation spectra from cosmic shear spectra at the first 4 redshift values 
use to integrate. This Figure is analogous to Figure 5, just extending to larger fc's. 
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Figure 10: As previous Figure, for 4 more redshift values. Also here, as in the Halofit case, 
discrepancies are smaller for larger z values. 



renormalize these spectra to obtain 

P R (k,z) =Af{k,z)P(k,z) 



(4.2) 



In Figure 11 some results of these operations are shown. They were obtained by ar- 
bitrarily choosing n s = 1. We however tested that the residual errors in the 
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Figure 11: Results of the inversion algorithm based on 12 point integration, if applied to "exact" 
spectra obtained from the simulation, by using the renormalization technique described in the text. 
The solid black curves (somewhere invisible) are the input fluctuation spectra P(k,z s ) (z s are the 
redshifts used by the Gaussian 12-point integration algorithm). The dotted curves are the spectra 
obtained by applying the inversion algorithm on the shear spectra. After renormalizing them, we 
obtain the solid red curves. For z >^ 0.7 the discrepancy keeps within 2-3% at all fc's. At lower 
z some spikes appear, indicating a failure of the SVD method at specific £ values. Below z ~ 0.5 
discrepancies mostly exceed 10 %. 



renormalized spectra are insensitive to the choice of n s in the interval 0.8 1.15. 
No tests were made outside this interval. 

At rather large z (> 0.7), the spectral details can be recovered, in this way. Therefore, 
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in principle, by re-exploiting the same information needed to create the window functions, 
i.e. by assigning Q m and h, we can recover parameters as as, H&, ^c (m.s.a. of fluctuations, 
baryon and CDM density parameters, respectively), as well as many baryon physics effects 
affecting large-£ shear spectra. At lower z, a precise spectral recovery is more difficult. 

We wish to add that the procedure still needs to be tested for the case when 
the DM component, or a part of it, is not cold. In particular, massive neutrinos, 
even with a mass in the 10 1 eV range, would cause spectral distorsions and 
might complicate the recovery. 



5. Discussion 
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The main achievement of this work is the development of a procedure to invert the equation 
yielding shear spectra Cij(£) from fluctuation spectra P(k,z), and its performance tested 
with spectra obtained either from Halofit or hydrodynamical simulations. 

As is known, lensing data yield the shear spectra Cij(l), while other observables yield 
the fluctuation spectrum P(k, z). Converting Cij(£) into P(k, z) would then be a valuable 
aid to compare the whole dataset with cosmological models. 

The stability of the inversion was tested here by applying the inversion matrix derived 
from a 12-point Gaussian integration to the "exact" Cij{£) spectra. The discrepancies 
between the latter Cij{£) and those obtainable through a 12-point Gaussian integration 
were as high as ~ 10%, in the worst points. It does not come as a surprise, then, that 
a brute application of the inverse "Gaussian matrix" does not yield satisfactory results. 
A renormalization procedure was however introduced, allowing us to recover fluctuation 
spectra for z >~ 0.7 with errors ~ 2-3% and down to z ~ 0.5 with errors ~ 10%. Lower 
redshift spectra are more difficult to recover. 

Such renormalization, requires that a few model parameters, defining the space-time 
"geometry" , which are assigned a priori. Besides recovering these parameters, the inversion 
then aims at constraining the "dynamical" parameters, e.g. as, Qf,, and, in addition, the 
effects of baryon physics. 

The restriction to 5 redshift bins and ~ 12 Gaussian points will be lifted, when Euclid 
data are available. With 7 (10) redshift bins, up to 28 (55) points are usable in Gaussian 
integration. In turn, the discrepancy between Gaussian and "exact" spectra scales as 
(n!) 2 /(2n)! ~ 2~ 2n [§5|. With 20 integration points, it is then 0(1O" 6 ) at most. Then this 
will no longer be the problem: applying SVD inverted matrix to "exact" results or, as 
matters most, to physical data should be satisfactory. 

The problem to be solved will then concern observational errors. We expect them to 
yield a sort of random noise. If so, SVD techniques provide the best method to deal with 
noisy matrix inversion. They work as better as redundancy increases. The point is whether 
a system of 55 equations dealing with 20-25 unknowns converges to the physical result, 
when data have a random noise 0(1%) or greater. 

In our opinion, perspectives are promising; the parameters beyond geometry could be 
recovered even in the presence of a systematic noise 10 times greater than the expected 
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noise in data. The main points to test, within this context, are whether the technique 
recovers spectra down to z ~ 0, and how data noise affects parameter precision. 

In the literature one finds cases when SVD techniques converge, but to spurious results. 
This is more and more unlikely as the system redundancy increases (see, e.g., [18|). We 



also recall that, in our case, the technique is to be applied to each i value separately. It is 
then even more unlikely that convergence to the same misleading result separately occurs 
at any £. We plan to devote further work to these issues. 

Admittedly, however, Euclid is a somewhat distant perspective. Using earlier available 
data, one could however try to achieve a better inversion, if galaxies could be shared among 
~ 7 bins. With the number of galaxies in datasets soon available, this might become 
possible if a substantial improvement in their redshift estimates is achieved, so that the 
bin dependence of Cy spectra stands out above shot noise. Euclid measures aim to reduce 
cr z /(l + z) from ~ 0.05 (the value used here) to ~ 0.03 . However, with fewer galaxies one 
could hope to do better, perhaps by measuring a spectroscopic redshift for a substantial 
fraction of the samples. 

Let us finally suggest that the stability problem might can be turned into an advantage. 
It is quite likely, in fact, that the inversion algorithm approaches regular results only when 
the choice of "geometrical" cosmological parameters is correct, yielding otherwise distorted 
results, e.g. an awkward k dependence. In association with the use of ~ 20 integration 
points, this could become a direct test to recover also background model parameters. 

The route opened by the inversion procedure discussed in this paper seems therefore 
promising, and might allow us a more complete exploitation of data coming from advanced 
cosmic shear measurements. 
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